Rotating thin-disk galaxies through the eyes of 
Newton 



James Q. Feng and C. F. Gallo 

CN Superconix Inc., 2440 Lisbon Avenue, Lake Elmo, MN 55042, USA 



Q* J E-mail: info@superconix.com 

< 



Abstract. By numerically solving the mass distribution in a rotating disk based 
on Newton's laws of motion and gravitation, we demonstrate that the observed flat 
rotation curves for most spiral galaxies correspond to exponentially decreasing mass 
density from galactic center for the most of the part except within the central core 
and near periphery edge. Hence, we believe the galaxies described with our model are 
consistent with that seen through the eyes of Newton. Although Newton's laws and 
C^ Kepler's laws seem to yield the same results when they are applied to the planets in 

ft the solar system, they are shown to lead to quite different results when describing the 

stellar dynamics in disk galaxies. This is because that Keplerian dynamics may be 

.£3 equivalent to Newtonian dynamics for only special circumstances, but not generally 

~2 for all the cases. Thus, the conclusions drawn from calculations based on Keplerian 

i— i dynamics are often likely to be erroneous when used to describe rotating disk galaxies. 
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1. Introduction 

A galaxy is a stellar system consisting of a massive gravitationally bound assembly of 
stars, an interstellar medium of gas and cosmic dust, etc. Observations have shown 
that many (mature spiral) galaxies share a common structure with the visible matter 
distributed in a flat thin disk (as in figure [I]), rotating about their center of mass in 
nearly circular orbits [1] . Apparently, this typical behavior of galaxies is similar to that 
of our solar system with planets orbiting around the Sun in a flat planetary plane. 




Figure 1. Spiral galaxies: NGC 3198 and NGC 4594-also known as M 104 sombrero 
galaxy. 

For planets orbiting around the Sun, Kepler's laws of planetary motion (obtained 
empirically) can provide accurate description. Yet it was Isaac Newton, in his 
"Philosophiae Naturalis Principia M athematica" , who used mathematical expressions 
to show that Kepler's laws are consequences of Newton's laws of motion and universal 
law of gravitation [2] . In addition, Newton found that Kepler's laws were only part of 
the story of how objects move in response to gravity. With his laws being discovered 
by analyzing the orbits of planets around the Sun, Kepler had no reason to believe that 
his laws would apply to other cases such as moons orbiting planets or comets orbiting 
the Sun. But Newton was able to derive more general rules that can explain the motion 
of objects throughout the universe, by analyzing his equations of gravity and motion. 
Therefore, Newton's laws of motion and gravity have become a crucial part of the 
foundation of modern astronomy, whereas Kepler's laws can become misleading if not 
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applied correctly with sufficient care to cases other than planets orbiting the Sun. For 
example, Kepler's third law, stating that more distant objects rotate around the center 
at slower average speeds, cannot describe the typically observed flat orbital velocity that 
remains invariant for most part of a galaxy outside its central core pH HI [6] . 

The fundamental difference between a galaxy and the solar system is that the mass 
is apparently distributed across the entire galaxy whereas the solar system has its mass 
concentrated at the center in the Sun. Each planet in the solar system can be quite 
reasonably treated as a point mass moving in a spherically symmetric gravitational field 
stemming from a large central point mass. The spherical symmetry of gravitational 
field in the solar system greatly simplifies the mathematical analysis, because the 
gravitational potential at any position is basically determined by the distance from 
the center and the mass of the Sun, with contributions from other planets negligible. 
Actually, the treatment similar to that for the solar system can be directly extended to 
the situation of distributed mass system as long as it retains the spherical symmetry, 
except that now the equivalent mass at the center is not a constant but equals to the mass 
enclosed by the concentric spherical surface through that point of interest. According to 
Newton's first and second theorems, if certain amount of mass is uniformly distributed 
in a spherical shell, this shell exerts no net (gravitational) force on any mass at any point 
inside it but attracts any mass outside it as if its mass is concentrated at the center 
p. Unfortunately, the thin-disk galaxies do not possess such a simplification-enabling 
spherical symmetry. So, much more sophisticated mathematical treatments are needed 
to correctly apply Newton's laws to the thin-disk galaxies. 

Here in this work, we demonstrate an effective numerical method for computing 
either mass density distribution for a given orbital velocity profile or vice versa by 
solving the governing equations based on Newton's laws for an axisymmetric thin-disk 
galaxy of finite size. We also quantitatively illustrate the possible misleading results due 
to incorrect application of Kepler's laws to the same thin-disk galaxy. In other words, 
we show that the observed behavior of disk galaxies can be described and explained by 
application of Newton's laws, but not by Kepler's laws that may only be regarded as a 
special case derived from Newton's laws for spherically symmetric gravitational field. 

2. Governing equations 

Because much of the mass of a galaxy resides in stars, we can in principle compute 
gravitational field of a large collection of stars by adding the point-mass fields of all 
the stars together for any spot of interest. For convenience of mathematical treatment, 
however, here we represent a galaxy by a continuum of axisymmetrically distributed 
mass in a circular disk of radius R g as shown in figure [2j Thus, in the present model we 
consider continuous mass density distribution instead of discrete mass points scattered 
throughout the disk. This kind of approximation is typically valid when the mass 
density in stars is viewed on a scale that is small compare to the size of the galaxy, 
but large compared to the mean distance between stars [1J. Physically, the stars in a 
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galaxy must rotate about the galactic center to maintain the disk shape. Without the 
centrifugal effect due to rotation, the stars would collapse into the galactic center as a 
consequence of the gravitational field. It is also reasonable to assume the galaxy is in an 
approximately steady state with the gravitational force and centrifugal force balancing 
each other, in view of the fact that most disk stars have completed a large number of 
revolutions [I]. 

Let's consider the force density on a test mass at (r, 9 = 0) generated by the 
gravitational attraction due to the summation (or integration) of a distribution of mass 
density p(r) at position described by the variables of integration (f, 9). Here the distance 
between (r, 9) and (f, 9) is (f 2 + r 2 — 2fr cos 9) x l 2 and the vector projection between the 
two points is (rcosO — r). Thus in a steady state, the mechanical balance between the 
gravitational force (due to the summation of mass in a series of concentric rings) and 
centrifugal force at every test point (r, 9 = 0) on the disk, according to Newton's laws 
of motion and gravitation, can be written as an integration equation (in terms of force 
per unit mass) 

"- 71 (rcos§-r)d§ 1 ,^, AJ . A V{r) 2 n 

— — ; — - — . p(r)hrdr + A w = 0, (1) 

(f 2 + r 2 — 2fr cos#) 3 / 2 r 

where all the variables are made dimensionless by measuring lengths (e.g., r, f, h) in 
units of the outermost galactic radius R g , disk mass density (p) in units of M g jR? g with 
M g denoting the total galactic mass, and velocities [V(r)] in units of the characteristic 
galactic rotational velocity V (as usually defined according to the problem of interest). 
The disk thickness h is assumed to be constant and small in comparison with the 
galactic radius R g . The results are expected to be insensitive to the exact value of 
this ratio as long as it is small. There is no difference in terms of physical meaning 
between the notations (r,9) and (f, 9); but mathematically the former denotes the 
independent variables in the integral equation (for 9 = 0) whereas the latter the variables 
of integration. The gravitational force represented as the summation of a series of 
concentric rings is described by the first term (double integral) while the centrifugal 
forces are described by the second term in ([T]). 

Our process of nondimensionalization of the force-balance equation yields a 
dimensionless parameter, which we call the "galactic rotation number" A, as given 
by 

A = V R 9 (2 ) 

A ~ M g G' [2} 

where G (= 6.67 x 1CT 11 [m 3 /(kg s 2 )]) denotes the gravitational constant, R g is the 
outermost galactic radius, and Vq is the characteristic velocity (which may be equated 
here to the maximum asymptotic rotational velocity of a disk galaxy). This galactic 
rotation number A simply displays the ratio of centrifugal forces to gravitational forces. 
For typical galactic values of R g , Vq, and M g we obtain A ~ 1.6 as will be shown in 
detail later. 
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Figure 2. Definition sketch of the thin-disk model considered in the present work. The 
mass is assumed to distribute axisymmetrically in the circular disk of uniform thickness 
h with a variable density as a function of radial coordinate r (but independent of the 
polar angle 0). 



When solving for the mass density p(r) with V(r) given, we need to impose an 
overall constraint such that the total mass of the galaxy M g is constant, namely, 

2n / p{f)hfdf = 1. (3) 

Jo 

This constraint due to the conservation of mass can actually be used to determine the 

value of galatic rotation number A. 

It is known that the integral with respect to 9 in (hj) can be written as [I] (pp. 

72-73) 



2- 



(rcos$ — r)d6 



E(m) 



K(m) 



r(f — r) r(f + r) 



(4) 



'o (r 2 + r 2 — 2fr cos#) 3 / 2 

where K(m) and E{m) denote the complete elliptic integrals of the first kind and second 
kind, with 

Anrnr 



m = 



Thus, ([I]) becomes 



(f + r)' 



E(m) 






K{m) 

f + r 



p(r)hrdr + -AV(rf = . 



(5) 



(6) 
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Equations ^ and (|3| can be used to determine the mass density distribution p{r) in 
the disk, the galactic rotation number A, and subsequently the total galactic mass M g , 
all from measured values of V(r), R c , R g and Vo- Seemingly complicated as it might 
be, these equations for p(r) actually constitute a linear mathematical problem that 
guarantees uniqueness of solutions. Conversely, these equations can also be used to 
determine the orbital velocity V(r) if the mass density distribution p{r) is known. This 
is a well-defined mathematical problem completely deducible from the available input 
data. 

Because our governing equations are derived according to Newton's laws, they must 
be applicable to the solar system. As an example for p{r) = 5(r)/(jrr) (namely, the 
Dirac delta function in two-dimensional polar coordinates), Q or ^ would yield 

--+AV{r) 2 = 0, (7) 

which is exactly the familiar formula for so-called Keplerian velocity based on Kepler's 
third law of planetary motion 



where Vq V(r) is the dimensional orbital velocity, M g is basically the mass of the Sun, 
and R g r the dimensional distance from the Sun. 

For a galaxy with mass distribution that is not spherically symmetric, simple closed- 
form analytical solutions may not be tractable. Yet accurate numerical solutions can 
be computed with appropriately implemented computational techniques as detailed in 
Appendix A. What it amounts to is nothing more than solving a linear algebra matrix 
problem using a well-established matrix solver. 

3. Mass distribution determined from given rotation curve 

The measurements of galactic rotational velocity profiles (also known as "rotation 
curves" ) of mature spiral galaxies reveal that the rotation velocity typically rises linearly 
from the galactic center (as if the local mass were in rigid body rotation) in a relatively 
small core, and then with the slope decreasing continuously in a narrow transition zone 
it reaches an approximately constant (flat) velocity extending to the galactic periphery 
[31 IH E] • These essential features may be mathematically idealized adf] 

V{r) = 1 - e- r/Rc , (9) 

where V(r) denotes the (dimensionless) orbital velocity (as measured in units of the 
maximum velocity in the flat part that may be regarded as the characteristic velocity 
of galactic rotation Vo), and r the radial coordinate from the galactic center (in units 

| Because the equations are solved numerically, the form of V(r) can be almost arbitrary. The 
idealized form of V(r) presented here is just for convenience of illustrations rather than the limitation 
of mathematical solution techniques. 
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Figure 3. Nondimensionalized orbital velocity profiles V(r) according to the 
mathematically idealized formula fcty for R c = 0.015, 0.03, 0.05 and 0.1. 



of the outermost galactic disk radius R g ). The parameter R c can be used to describe 
various radii of the "cores" of different galaxies. Close to the galactic center, namely 
when r/R c is small, we have V(r) ~ r/R c , describing a linearly rising rotation velocity, 
by virtue of the Taylor expansion of e~ r ^ Rc . Figure |3j shows typical galactic rotation 
curves described by ^. 

With V(r) given by (J9J), the mass density distribution p{r) and the value of A 



can be determined by computing solution to (A. 3) and (A.5). To compute numerical 
solutions, the value of disk thickness h must be provided; we assume h = 0.01 (based 
on the measurements for the Milky Way galaxy). By using N = 1001 nonuniformly 
distributed nodes we found the obtained curves of p(r) become reasonably smooth and 
the values of galactic rotation number A are discretization- independent. 

Shown in figure [I] are the computed mass density distributions p(r) that satisfy the 
galactic rotation curves in figure |3j It is at the galactic center where mass attains the 
highest density. Away from the galactic center, the mass density decreases rapidly (with 
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Figure 4. The distributions of mass density p(r) computed based on Newtonian 
dynamics and given rotation curves for R c = 0.015, 0.03, 0.05, and 0.1, with 
A = 1.5714, 1.5733, 1.5777, and 1.5999 determined as part of the numerical solutions. 



a slope becoming steeper for a tigher galactic core indicated by smaller R c ). However, 
beyond R c , the mass density decreases rather gradually towards the galactic periphery 
until reaching the galactic edge where it takes a sharp drop. Noteworthy here is that 
the computed values of galactic rotation number A are within a small range around 1.6 
despite an order-of-magnitude variation of the galactic core radius R c . 

As apparent in figure [4j the computed log p decreases almost linearly with r except 
for r < 0.1 (within the central galactic core) and r > 0.9 (near the galactic edge). This 
general feature is quite similar to the measured brightness distributions (in typical spiral 
galaxies) that are commonly fitted in an exponential form with regions of central core 
and outer edge truncated. For the case of R c = 0.015 and A = 1.5714, a least-square fit 
of our computed logp versus rfor0.1<r<0.9 yields 

logp = 5.4179- 3.6802 r. (10) 
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This actually corresponds to an exponential function p = po e~ r l Rd with 

p = 225.4 and R d = 0.2717 , (11) 

which is known not to be able to generate the observed flat rotation curves [7J, Q]. Thus 
the much more rapid decrease of logp in the small intervals [0, 0.1) and (0.9, 1] must 
play important roles in compensating the deficiencies of the simple exponential mass 
density distribution for matching the commonly observed (flat) rotation curves. If we 
assume the mass density follows roughly the exponential distribution p = p e~ r / Rd , we 
would have the cumulative mass from the galactic center given by 

rr 

2nh p e- f/Rd rdr = 2Trhp [R 2 d -R d (r + R d )e- r/Rd }. (12) 

Jo 

At r = 1, this yields a value of 0.922 for p = 225.4 and Rd = 0.2717, indicating the 

exponential mass density distribution is likely to describe about 90% of the mass in a 

disk galaxy. Another ~ 10% of the galactic mass seems to reside in the galactic core 

which may not follow the same exponential distribution (cf. figure El). 

From the knowledge of Vq and R g from measured rotation curves, we can determine 

the value of M g based on computed value of the galactic rotation number A (cf. ([2])) as 



M,M. (13) 



Because our computed A varies little from 1.6 for various R c , (13) suggests a general 
relationship of M g oc V 2 R g as what Bosnia j5] found from evaluating mass versus size 
in a large number of observed disk galaxies. In view of the fact that the values of Vq are 
typically around 200 (km/s), we then believe that the disk size of galaxies R g must be 
finite in order to keep the total mass of a galaxy M g from becoming infinity. Suggesting 
finite disk size of galaxies does not necessarily mean that the mass density becomes 
zero outside the disk edge. We believe that the mass density beyond the galactic edge 
approaches the inter-galactic level of value and is roughly spherically symmetric, which 
leads to inconsequential gravitational effect on rotation dynamics in the disk part. 

If we take the rotation curve with R c = 0.015 in figure [3] as that of the Milky Way, 
we have the galactic rotation number A = 1.57[§J Then, from measured Milky Way 
values V Q = 2.2 x 10 5 (m/s) and R g = 5x 10 4 (light-years) = 4.73 x 10 20 (m), Q yields 

M g = 2.19 x 10 41 (kg) = 1.10 x 10 11 (solar-mass) . (14) 

This value is in very good agreement with the Milky Way star counts of 100 billion [8], 
further including additional dust, grains, lumps, gases and plasma in all galaxies. 

With the given values of M g and R g , we can estimate the computed 'radial scale' 
RdRg = 4.0 (kpc) = 1.234 x 10 20 (m) and the exponential disk central (surface) mass 



density pohM g jR g = 1340 (solar-mass / pc ) based on (10) for the Milky Way galaxy. 
(Here, 1 (pc) = 3.086 x 10 16 (m) and 1 (solar mass) = 1.99 x 10 30 (kg).) Compared 
with the results from fitting the brightness measurement data, e.g., the radial scale 

§ Here we use the case of R c = 0.015 only for the purpose of convenient illustration. The value of 
A = 1.57 (with a variation of 1%) is actually valid for a wide range of R c values as shown in figure El 
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of 2.5 (kpc) and exponential disk central brightness of 867 (solar-luminosity / pc 2 ) 
[9], our computational results indicate either generally dimmer stars (than the Sun) or 
considerable amounts of cold gas exist throughout the Milky Way with the total mass 
density decreasing at slower rate than that of the brightness (i.e., luminosity density). 
Another example is the galaxy NGC 3198, which has a (nearly idealized, often 
cited) rotation curve with V = 1.5 x 10 5 (m/s), R g = 30 (kpc) = 9.24 x 10 20 (m), and 
R c ~ 0.015. Again with A = 1.57, we obtain M g = 1.98 x 10 41 (kg) = 9.9 x 10 10 (solar- 
mass). As with the Milky Way, we can also predict the radial scale and exponential 
central mass density for NGC 3198 as 8.16 (kpc) and 250 (solar- mass / pc 2 ), respectively 



(based on p = 225.4 and Rd = 0.2717 from (10). Compared with the radial luminosity 
profile (which suggested an exponential disk with a radial scale of 2.63 (kpc) and central 
brightness 212 (solar-luminosity / pc 2 ) based on a total luminosity of 9.0 x 10 9 (solar- 
luminosity) [10] , our predicted mass density appears to decrease much more slowly with 
stars (or mass objects) generally dimmer than the Sun. 

Hence, the results computed here, as if they were obtained by Newton applying 
his laws of motion and gravitation to solve the governing equations (flj and (pi), seem 
to be reasonably consistent with the observational measurements. In other words, the 
rotating thin-disk galaxies through the eyes of Newton are nothing more than massive 
gravitationally bound assemblies of objects governed by his same laws for the planet's 
motion in the solar system albeit more sophiscated mathematical treatments are needed 
to obtain the correct description than those with the planets in the solar system. 

4. Results based on Keplerian dynamics 

In the literature, many authors (21 El [HI H2J EE2 EH] tend to simply apply Keplerian 
dynamics (which was derived from gravitational field generated by a spherically 
symmetric distribution of mass) when analyzing the rotation behavior of thin-disk 
galaxies. For example, Volders [H] demonstrated that spiral galaxy M33 does not 
spin as expected according to Keplerian dynamics-a result which was later extended to 
many other spiral galaxies [21 IH El E] • 

Strictly speaking, a Keplerian potential (due to a point mass M as that for the 
solar system) is expressed as 

$( r ) = (15) 

r 

For a distributed mass with spherical symmetry, the generalized form of Keplerian 

potential becomes 

GM(r) , s 

$ r = U 16 

r 
where M(r) denotes the amount of mass enclosed by the concentric spherical surface of 



radius r [lj. Although (16) comes from the assumption of spherical symmetry, it has 
often been used to determine the mass distribution and the total mass of a (disk) galaxy 
from a measured rotation curve, with M(r) denoting the mass interior to radius r from 
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the galactic center. For example, in the recent versions of textbooks by Bennett et al. 
[2], Sparke and Gallagher [8], and Keel [12], the value of M{r) (also denoted as M r or 
M(< r)) in a disk galaxy is simply determined from a known rotation curve V(r) by 

rV(r) 2 



M(r) 



G 



(17) 



which has basically the same mathematical form as pi). 



According to (17), we would have an equation based on Keplerian dynamics for 
force balance as 

0, (18) 



2n h / p{r)rdr — ArV{r)' 



instead of (pi). Here, the difference in mathematical forms between (18) and (ll| should 
be quite clear. But whether there are much of practical differences between the two 
may not be obvious. For example, authors like Sparke and Gallagher [8] and Keel [12] 
attempted to justify their usages of Keplerian formula for rotating thin disk galaxies by 
stating that the result due to Keplerian formula does not differ more than 20% from 
the actual result, without showing a quantitative comparison. Therefore, we would like 



to examine a few comparative examples to see whether Keplerian dynamics ( 18 ) can be 



practically used as a close approximation to Newtonian dynamics ([!]) for disk galaxies. 
For the orbital velocity V(r) given by ([9]), an analytical solution to (18) for p(r) 
can be obtained as 



p{r) 



A 

2^h 



(1 _ 2e"^ + e ^r/B^ + ±_ ( e -r/R c _ g-Sr/flc) 



2 

R r 



(19) 



It is not difficult to prove that p(r) — > as r — >■ 0, as in sharp contrast to that obtained 
according to Newton's laws shown in figure 111 For r >> R c (i.e., outside the galactic 



core), (19 describes p{r) oc 1/r, as expected for the part where V(r) is flat (constant) 



such that the first term in (18) becomes proportional to r. In fact, many astronomers 
often consider the flat rotation curves, i.e., V(r) = constant, to indicate that the mass 
of a disk galaxy should continue to increase linearly well beyond its bright central region 
because r V{r) 2 oc r 



To satisfy (|3]), the value of A in (19) is given by 

1 



A 



1 



When R c is small, e.g. 



. 2e-V«c + e -2/n c ■ 
R c = 0.015, we have e" 



(20) 



X I R ° -> 0, e.g., e" 1 / - 015 ~ 10~ 29 ; therefore, 
A m 1 according to (20). This is in sharp contrast to the result based on formulas for 
disk galaxies (which predicts A = 1.57). If A = 1 were used in (13), we would have the 
total mass in Milky Way equal to 2.2 x 10 11 (solar-mass) -much more than that given 



by (14) and the Milky Way star counts. 



As a comparison, the distribution of p(r) shown in figure for R c = 0.015 and 
that given by ( 19 ) with A = 1 are shown in figure [5J Now, the differences between the 



two are obvious: the Keplerian mass density behaves totally differently in the galactic 
core with p(r) — ¥ whereas that numerically obtained in § 3 has p(r) monotonically 
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Figure 5. The distribution of mass density p(r) computed based on Newtonian 
dynamics for R c — 0.015 (solid) and that derived based on Keplerian dynamics given 



by (19) with A = 1 and R c = 0.015 (dotted). 



decreasing with r from its maximum value at galactic center to periphery. Outside the 
galactic core, the Keplerian mass density shows much slower decay (oc 1/r) toward the 
galactic periphery whereas that obtained numerically in § 3 seems to be approximately 
exponential. Therefore, the mass density distribution determined based on Keplerian 



dynamics (18) cannot be a good approximation to the actual p(r) obtained numerically 
in § 3 with Newton's laws strictly applied. 

On the other hand, if the mass density p(r) is known as obtained in § 3 for the 
case of R c = 0.015 with A = 1.5714, we can compute the corresponding V(r) from (18). 



figure M shows that the V(r) from (18) clearly differs from that of ([9j, with rotation 
velocity decreasing with the radial distance r toward the galactic periphery as expected 



from Kepler's third law. Again, this suggests that (18) according to the Keplerian 



dynamics cannot be a good approximation to the actual galactic dynamics. 

As demonstrated in § 3, we believe that Newton's laws of motion and gravitation 
can adequately describe the dynamical behavior of (disk) galaxies, with appropriate 
mathematical treatments. Kepler's laws should not be regarded as the same as Newton's 
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Figure 6. The rotation curve of R c = 0.015 described by pi (solid) from Newtonian 
dynamics and that determined from ( 18 1 based on the Keplerian dynamics with the 



known p(r) for R c = 0.015 obtained in § 3 with A = 1.5714 (dotted). 

laws. Newton's laws can explain Kepler's laws for the planets in the solar system; but 
Kepler's laws cannot be extended to the galactic dynamics like Newton's laws, even in 
an approximation sense as shown in figure [5] and figure [6j 

5. Conclusions 

By strictly applying Newton's laws, the present computational model can yield mass 
distributions from the observed galactic rotation curves, as apparently quite consistent 
with the observed near exponential brightness distributions. Thus through the eyes 
of Newton, galaxies are nothing more than graviationally bound assemblies of massive 
objects that are governed by his same laws for the planet's motion in the solar system. 
Although Newton's laws and Kepler's laws seem to yield the same results when they are 
applied to the planets in the solar system, they can lead to quite different results when 
describing the stellar dynamics in disk galaxies (cf. figure [5] and figure pi). If not careful, 
simply extending Kepler's laws to the disk galaxies with subtensively distributed mass 
can suggest misleading conclusions. 
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As demonstrated in § 4, substituting the computed mass density distribution p{r) 



based on Newtonian dynamics into the Keplerian force-balance equation (18) would 
yield a rotation curve with orbital velocity decreasing toward galactic periphery (see 
figure [6]). This had led many authors to believe that the visible mass in a galaxy 
cannot explain the observed flat rotation curve [El El E]- Therefore, some authors 
have speculated that some kind of invisible matter called dark matter must exist in 
the galaxy [21 El El E]- Other authors believed modification of Newton's laws to 
be needed [IB] . The fundamental problem here is that astronomers tend to determine 
the 'visible' mass in a galaxy from the measured brightness based on an over-simplified 
mass-to-light ratio [15J. When the mass distribution so estimated did not generate 
the observed flat rotation curve, especially when the Keplerian dynamics was used for 
another over-simplification, it is often referred to as the 'galactic rotation problem' 
suggesting that there is a discrepancy between the observed rotation speed of matter 
in the disk and the predictions of Newtonian dynamics [21 El El El E]- But we 
believe that Newtonian dynamics can adequately explain the stellar dynamics in disk 
galaxies when applied correctly; neither the introduction of dark matter nor modification 
of Newtonian dynamics is needed for explaining the observed rotation curve [T71 [PS] . 
Hence, the rotating disk galaxies described with our model must be that seen through 
the eyes of Newton. 

Appendix A. Computational techniques 

Following a standard boundary element method [121 120], the governing equations ffity 
and (|3]) can be discretized by dividing the one-dimensional problem domain [0, 1] into a 
finite number of line segments called (linear) elements. Each element covers a subdomain 
confined by two end nodes, e.g., element n corresponds to the subdomain [r n ,r n +i], 
where r n and r n+1 are nodal values of r at nodes n and n + 1, respectively. On each 
element, which is mapped onto a unit line segment [0,1] in the ^-domain (i.e., the 
computational domain), p is expressed in terms of the linear basis functions as 

p(0 = p n (l - + p n+1 £ , < £ < 1 , (A.l) 

where p n and p n+ \ are nodal values of p at nodes n and n + 1, respectively. Similarly, 
the radial coordinate r (as well as f) on each element is also expressed in terms of the 
linear basis functions by so-called isoparameteric mapping: 

r(0=r n (l-0+r„+i£, 0<£<1. (A.2) 

If V(r) is given (e.g., from measurements), the N nodal values of p n = p(r n ) can be 
determined by solving N independent residual equations over N — 1 element obtained 
from the collocation procedure, i.e., 



N-l 

E 

71=1 



E{rrii) K(rrii 



r(£)-n r(£) + ri 



dv ! 

p(0hm^+-AV(n) 2 = , i = 1, 2, ..., N , (A.3) 
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with 



mi(£) = 



4f(£)r; 



[f(0+r,] 2 ' 

where p(£) = p n (l - £) + Pn+i£ and r(£) = r„(l - 
be solved by the addition of the constraint equation 

N-l 



— i' 1 dr 

n=l 



d£ 



(A.4) 
r n +i£. The value of A can also 

(A.5) 



Thus, we have iV + 1 independent equations for determining AT + 1 unknowns. The 



mathematical problem is now well-posed. The set of linear equations (A. 3) and (A.5) 



for N + 1 unknowns, i.e., N nodal values of p n and A, can be written in a matrix form 

as 

Jx = -R, (A.6) 

where R is the residual vector consisting of N + 1 components given by the left side of 



(A. 3) and (A.5), x is the unknown vector of N nodal values of p n and A, and J is the 
Jacobian matrix of sensitivities of the residual R to the unknowns x, i.e., J^ = dRi/dxj. 



The matrix equation (A.6) is actually derived based on Newton's method (also known 
as the Newton-Raphson method) for iteratively finding roots of a set of multi-variable 
nonlinear functions. For a set of linear functions, as in the present case, a single iteration 
is enough for obtaining the solution. 

The complete elliptic integrals of the first kind and second kind can be numerically 
computed with the formulas [21] 

4 4 



K(m) = 2, Q>iTn\ — log(mi) 2_, bim\ 



1=0 



1=0 



and 



where 



E{m) = 1 + 2_. c i m [ ~ l°g( m i) /_] d 



I'm' 



i • 



i=i 



i=i 



mi 



m 






f + r 



(A.7) 



(Ai 



(A.9) 



Thus, the terms associated with Kirrii) and Eirrii) in (A.3) become singular when f — y r, 
on the elements with r« as one of their end points. 

The logarithmic singularity is treated by converting the singular one-dimensional 
integrals into non-singular two-dimensional integrals by virtue of the identities: 

/ 1 /(e)iogede=-/o 1 /o 1 /(^)^de 

j l /(o io g (i - ode = - £ £ /(i - fr)*7de ' 

where /(£) denotes a well-behaving (non-singular) function of £ on < £ < 1, which 
can be derived by considering integration over a triangular area in a two-dimensional 



(A.10) 
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fit) 



X 



d£ ) dx 



1 x J 



1 



m) 



X 



dx I d£ 



f(r)x) xdrj I dx 



1 rl 



JO 



f(r)£)dr)d£. 



and 



/(O io g (i - o ^ 



1-!/ 



/(fl 



o vi 

1 /•! 



d£ dj/ 



i 



i-« 



/(O 



dy ^ 



, , f(i-vy)ydr)) d y 

o y \Jo J 



JO 



f(l-r)£) drift. 



But a more serious non-integrable singularity l/(f — r$) exists due to the term 
E{rrii)/(r — ri) in (A.3) as r — > r^. The l/(f — r») type of singularity is treated by taking 
the Cauchy principal value to obtain meaningful evaluation [22]. In view of the fact 
that each T{ is considered to be shared by two adjacent elements covering the intervals 
b"i-i, r i] an d [ r «; r i+i]i the Cauchy principal value of the integral over these two elements 
is given by 



lim 

e->0 



"* e p(f)fdf f n+1 p(f)fdf 



L-zn-x 



r — Tj 



r — Ti 



(A.11) 



In terms of elemental £, (A. 11) is equivalent to 

'W(n-r«-i) [^(1-0 + ^[r w (l - + r^de 
'o 



lim 

e-*0 



1-^ 



e/(r i+1 -n) 



\Pi(i - o + pmeiM - o + n+ieiji ] 



(A.12) 



Performing integration by parts on (|A.12|) yields 

V,-,i - r. 



Pi Ti log 

+ 



i+l 



1 d{^_x(l - + PiUn^l - + r£\} 



n - n-x j \j Q dt, 

1 d{\p i (l-0+Pi+i^][n(l-0+ri + ^]} 



io g (i - m 



log^ 



'0 ^ 

where all the terms associated with loge cancel out each other, the terms with eloge 
become zero at the limit of e — > 0. The first term becomes nonzero when the mesh 
notdes are not uniformly distributed (namely, the adjacent elements are not of the same 
segment size). 

At the galaxy center n — 0, 

-n+i p(f)rdr rn+1 



r — Vi 



p(f)df . 



Thus, the l/(f — r t ) type of singularity disappears naturally. 



Rotating thin-disk galaxies through the eyes of Newton 17 

When r% — 1 (i.e., i = N), it is the end node of the domain. We can use 
a numerically relaxing boundary condition by imagining another element extending 
beyond the domain boundary covering an interval [r», ri+i], because it is needed for 
the treatment with Cauchy principal value. In doing so we can also have r^x — r^ = 
Ti — Ti_i such that log[(r i+1 — rj)/(fj — rj_i)] becomes zero, to simplify the numerical 
implementation. Moreover, it is reasonable to assume that pi + \ = because it is located 
outside the disk edge where the extremely low intergalatic mass density is expected to 
have inconsequential gravitational effect. With sufficiently fine local discretization, this 
extra element can be considered to cover a diminshing physical space such that its 
existence becomes numerically inconsequential. Thus, at r, = 1 we have 

11 d{[ Pi (i -0+ A+i£]h(i - + n+i£]} 



{pi+i - Pi) I r(0 log£d£ + (r i+1 - n) I p(0 log^ = p t 



Ti ~ -[Ti -Ti-\) 



Now that only logarithmic singularities are left, (A. 10) can be used to eliminate all 
singularities in integral computations. 

In addition, to avoid cusps in mass density at the galactic center, continuity of the 
derivative of p at the galaxy center r = is applied when solving for p with given V(r). 
This boundary condition is imposed at the first node % = 1 to require dp/dr = at 
r = 0, which becomes 

p(r 1 )=p(r 2 ) (A.13) 

in discretized form. 

Noteworthy here is that the (removable) singularities in the kernels of the integral 
equation (pi), when properly handled, lead to a diagonally dominant Jacobian matrix 



in (A.6) with bounded condition number. This fact makes the matrix equation (A.6) 
quite robust for almost any straightforward matrix solvers. In the present work, we 
simply used the available code for Gauss elimination [23]. To check the correctness of 
our computational code implementation, we substituted an exponential mass density 
distribution (e.g., p(f)h = e~ 5r ) into ^ and compared the computed orbital velocity 
V{r) with the well-known analytical formula of Freeman [7]. The result showed excellent 
agreement. Moreover, we could also obtain constant orbital velocity V(r) = 1 by 
integrating the density of Mestel's disk (21] p{r) = A[l — (2/7r) sin _1 (r)]/(27r/w) through 

The numerical method presented here can be applicable to rotation curves of 
arbitrary forms and does not require assumptions about the rotation curve beyond the 
radial coordinate where the orbital velocity is no longer measurable as in Ref. [1] when 
using the formula of Toomre [25] that contains an integral extending to infinity. Not 
only is it convenient for considering galaxies of finite disk sizes, it can also become an 
effective tool for deducing the mass distribution in a thin-disk galaxy from the measured 
rotation curve based on Newtonian dynamics. 
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